home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / raytrace / renaisnc / lib.lha / Lib / 4D.c next >
Encoding:
C/C++ Source or Header  |  1992-01-04  |  7.9 KB  |  451 lines
  1. /*
  2.   File: 4D.c
  3.   Author: K.R. Sloan, Jr.
  4.   Last Modified: 11 September 1985
  5.                  22 May 1986 by Tony DeRose
  6.                  17 September 1986 by Jamie Painter
  7.          04 November 1986 by Karen Foltz
  8.   Purpose: 4D homogeneous coordinate manipulation package
  9.  */
  10.  
  11. #include <stdio.h>
  12.  
  13. #include "Math.h"
  14. #include "4D.h"
  15.  
  16. #include "fast.c"
  17.  
  18. typedef short bool;
  19. extern bool Kdebug;
  20.  
  21.  
  22. /* Exports */
  23. PointType4D Origin = {0.0,0.0,0.0,1.0};
  24. double epsilon = 0.00000000001;
  25.  
  26. PointType4D CreatePoint4D(x,y,z,w)
  27. double x,y,z,w;
  28. {
  29.   PointType4D NewPoint;
  30.   NewPoint.P[0] = x;
  31.   NewPoint.P[1] = y;
  32.   NewPoint.P[2] = z;
  33.   NewPoint.P[3] = w;
  34.   return NewPoint;
  35. }
  36.  
  37. /*
  38. **  Return the point 1 t-th along the way
  39. **  from P0 to P1.
  40. */
  41. extern PointType4D Ratio4D(P0,P1,t)
  42. PointType4D P0,P1;
  43. double t;
  44. {
  45.   PointType4D P;
  46.  
  47.   P.P[X] = P0.P[X] + t * (P1.P[X] - P0.P[X]);
  48.   P.P[Y] = P0.P[Y] + t * (P1.P[Y] - P0.P[Y]);
  49.   P.P[Z] = P0.P[Z] + t * (P1.P[Z] - P0.P[Z]);
  50.   P.P[W] = 1.0;
  51.  
  52.   return P;
  53. }
  54.  
  55.  
  56. extern PointType4D PxT4D(P,T)
  57. PointType4D P;
  58. TransformType4D T;
  59.  {
  60.   PointType4D Q;
  61.  
  62.   PxT4Dm(Q,P,T);
  63.  
  64.   return Q;
  65.  }
  66.  
  67.  
  68. extern TransformType4D Transpose4D(T)
  69. TransformType4D T;
  70. {
  71.   TransformType4D Tt;
  72.   int i,j;
  73.  
  74.   for (i=0; i<4; i++) {
  75.     for (j=0; j<4; j++) {
  76.       Tt.T[i][j] = T.T[j][i];
  77.     }
  78.   }
  79.   return Tt;
  80. }
  81.  
  82.  
  83.  
  84.  
  85. extern TransformType4D TxT4D(A, B)
  86.  TransformType4D A, B;
  87.  {
  88.   TransformType4D C;
  89.  
  90.   TxT4Dm(C,A,B);
  91.   return C;
  92.  }
  93.  
  94. extern TransformType4D Identity4D()
  95.  {
  96.   static TransformType4D I = { 1., 0., 0., 0.,
  97.                                0., 1., 0., 0.,
  98.                                0., 0., 1., 0.,
  99.                                0., 0., 0., 1.};
  100.   return I; 
  101.  }
  102.  
  103. extern TransformType4D Translate4D(V)
  104. VectorType4D V;
  105. {
  106.   TransformType4D T;
  107.   
  108.   T = Identity4D();
  109.   T.T[0][3] = V.V[0];
  110.   T.T[1][3] = V.V[1];
  111.   T.T[2][3] = V.V[2];
  112.   return T;
  113.  }
  114.  
  115. TransformType4D UniformScale4D(s)
  116. double s;
  117. {
  118.   return Scale4D(CreatePoint4D(s,s,s,1.0));
  119. }
  120.  
  121. extern TransformType4D Scale4D(P)
  122.  PointType4D P;
  123.  {
  124.   TransformType4D T;
  125.  
  126.   T = Identity4D();
  127.   T.T[0][0] = P.P[0];
  128.   T.T[1][1] = P.P[1];
  129.   T.T[2][2] = P.P[2];
  130.   T.T[3][3] = P.P[3];
  131.   return T;
  132.  }
  133.  
  134. extern TransformType4D Rotate4D(P0, P1, Radians)
  135.  PointType4D P0, P1;
  136.  double Radians;
  137.  {
  138.   TransformType4D T, R1, R2, R3, Rotate, Tmp;
  139.   double w, a, b, c, l, v;
  140.   
  141.   w = P1.P[3]; a  = P1.P[0]/w; b  = P1.P[1]/w; c  = P1.P[2]/w;
  142.   w = P0.P[3]; a -= P0.P[0]/w; b -= P0.P[1]/w; c -= P0.P[2]/w;    
  143.  
  144.   l = sqrt(a*a + b*b + c*c);
  145.   if (l < epsilon)
  146.    {
  147.     fprintf(stderr,"Rotate4D: |P1-P0| is too small - can't rotate\n");
  148.     Tmp = Identity4D();
  149.     return Tmp;
  150.    }
  151.   else
  152.    {
  153.     a = a/l; b = b/l; c = c/l;
  154.    }
  155.   
  156.   v = sqrt(b*b + c*c);
  157.  
  158.   P0.P[3] = -P0.P[3];  /* reverse direction by reversing w */
  159.   T = Translate4D(Pdiff(P0,Origin)); 
  160.  
  161.   R1 = Identity4D(); R2 = R1; R3 = R1;
  162.  
  163.   if (v > epsilon)
  164.    {
  165.                           R1.T[1][1] = c/v;         R1.T[2][1] = b/v;
  166.                           R1.T[1][2] = -R1.T[2][1]; R1.T[2][2] = R1.T[1][1];
  167.  
  168.    }
  169.  
  170.   if ((v > epsilon) || (a > epsilon) || (a < (-epsilon)))
  171.    {
  172.     R2.T[0][0] = v;                                 R2.T[2][0] = a;
  173.     R2.T[0][2] = -R2.T[2][0];                       R2.T[2][2] = v;
  174.    }
  175.  
  176.   R3.T[0][0] = cos(Radians); R3.T[1][0] = -sin(Radians);
  177.   R3.T[0][1] = -R3.T[1][0];  R3.T[1][1] = R3.T[0][0];
  178.  
  179.   TxT4Dm(Rotate,T,R1);
  180.   TxT4Dm(Tmp,Rotate,R2);
  181.   TxT4Dm(Rotate,Tmp,R3);
  182.  
  183.                                                     R2.T[2][0] = -R2.T[2][0];
  184.   R2.T[0][2] = -R2.T[0][2];
  185.  
  186.                                                     R1.T[2][1] = -R1.T[2][1];
  187.                            R1.T[1][2] = -R1.T[1][2];
  188.  
  189.   T.T[0][3] = -T.T[0][3];   T.T[1][3] = -T.T[1][3];   T.T[2][3] = -T.T[2][3]; 
  190.  
  191.   TxT4Dm(Tmp,Rotate,R2);
  192.   TxT4Dm(Rotate,Tmp,R1);
  193.   TxT4Dm(Tmp,Rotate,T);
  194.   return Tmp;
  195.  }
  196.  
  197. extern TransformType4D Perspective4D(f)
  198.  double f;
  199.  {
  200.   TransformType4D P;
  201.   
  202.   P = Identity4D();
  203.   if ((f < epsilon) && (f > (-epsilon)))
  204.    {
  205.     fprintf(stderr,"Perspective4D: f (%f) is too small.\n",f);
  206.     if (f < 0.0) f = (-epsilon);
  207.     else         f =   epsilon;
  208.    }
  209.   P.T[2][2] = 1.0/f;    P.T[3][2] = 1.0/f;
  210.   P.T[2][3] = -1.0;     P.T[3][3] = 0.0;
  211.   return P;
  212.  }
  213.  
  214. void PrintTransform4D(f,T)
  215. FILE *f;
  216. TransformType4D T;
  217. {
  218.   int i,j;
  219.  
  220.   for (i=0; i<4; i++) {
  221.     fprintf(f,"[");
  222.     for (j=0; j<4; j++) {
  223.       fprintf(f, " %lg ", T.T[j][i]);
  224.     }
  225.     fprintf(f,"]\n");
  226.   }
  227. }
  228.  
  229. void PrintPoint4D(f,P)
  230. FILE *f;
  231. PointType4D P;
  232. {
  233.   fprintf(f,"[ %lg %lg %lg %lg ]", P.P[0], P.P[1], P.P[2], P.P[3]);
  234. }
  235.  
  236. /*
  237. ** Vectors are represented as 4-tuples with the last component
  238. ** being 0.  This representation allows them to be correctly transformed
  239. ** by homogeneous matrices.
  240. */
  241. VectorType4D CreateVector4D(a,b,c)
  242. double a,b,c;
  243. {
  244.   VectorType4D NewVector;
  245.  
  246.   NewVector.V[0] = a;
  247.   NewVector.V[1] = b;
  248.   NewVector.V[2] = c;
  249.   NewVector.V[3] = 0.0;
  250.   return NewVector;
  251. }
  252.  
  253. void PrintVector4D(f,v)
  254. FILE *f;
  255. VectorType4D v;
  256. {
  257.   fprintf(f,"<%lg,%lg,%lg>\n", v.V[0],v.V[1],v.V[2]);
  258. }
  259.  
  260. extern VectorType4D VxT4D(V,T)
  261. VectorType4D V;
  262. TransformType4D T;
  263.  {
  264.   VectorType4D Q;
  265.   int i,k;
  266.   double Qi;
  267.  
  268.   for (i = 0; i < 4; i++)
  269.    {
  270.     Qi = 0.0;
  271.     for (k = 0; k < 4; k++)
  272.      Qi += V.V[k] * T.T[i][k];
  273.     Q.V[i] = Qi;
  274.    }
  275.   return Q;
  276.  }
  277.  
  278. VectorType4D NormalxT4D(N,Tinv)
  279. VectorType4D N;
  280. TransformType4D Tinv;
  281.  {
  282.   VectorType4D Np;
  283.   int i,k;
  284.   double Npi;
  285.  
  286.   for (i = 0; i < 3; i++)
  287.    {
  288.     Npi = 0.0;
  289.     for (k = 0; k < 3; k++)
  290.      Npi += N.V[k] * Tinv.T[k][i];
  291.     Np.V[i] = Npi;
  292.    }
  293.   Np.V[3] = 0.0;
  294.   return Np;
  295.  }
  296.  
  297. VectorType4D Cross4D( v1, v2)
  298. VectorType4D v1, v2;
  299. {
  300.   VectorType4D v;
  301.  
  302.   v.V[0] = v1.V[1] * v2.V[2] - v1.V[2] * v2.V[1];
  303.   v.V[1] = -(v1.V[0] * v2.V[2] - v1.V[2] * v2.V[0]);
  304.   v.V[2] = v1.V[0] * v2.V[1] - v1.V[1] * v2.V[0];
  305.   v.V[3] = 0.0;
  306.  
  307.   return v;
  308. }
  309.  
  310. double Dot4D(v1,v2)
  311. VectorType4D v1, v2;
  312. {
  313.   return v1.V[0] * v2.V[0] + v1.V[1] * v2.V[1] + v1.V[2] * v2.V[2];
  314. }
  315.  
  316. double Vmag( v)
  317. VectorType4D v;
  318. {
  319.   return sqrt( v.V[0]*v.V[0] + v.V[1]*v.V[1] + v.V[2]*v.V[2]);
  320. }
  321.  
  322. VectorType4D Normalize( v)
  323. VectorType4D v;
  324. {
  325.   VectorType4D vout;
  326.   double l = sqrt( v.V[0]*v.V[0] + v.V[1]*v.V[1] + v.V[2]*v.V[2]);
  327.  
  328.   if (l < epsilon) {
  329.     fprintf(stderr, "Can't normalize a zero vector.\n");
  330.     /*exit(-1);*/
  331.     vout = v;
  332.   } else {
  333.     vout.V[0] = v.V[0]/l;
  334.     vout.V[1] = v.V[1]/l;
  335.     vout.V[2] = v.V[2]/l;
  336.     vout.V[3] = 0.0;
  337.   }
  338.  
  339.   return vout;
  340. }
  341.  
  342. /*
  343. ** Set W to 1
  344.  */
  345. #ifdef NOT_INLINE
  346. extern void NormalizePoint(P1)
  347. PointType4D *P1;
  348. {
  349.   register double *p = P1->P;
  350.   register double  w = p[W];
  351.  
  352.   if (w != 1.0) {
  353.     /* We are assuming that the order in P is:  X, Y, Z, W */
  354.     *p++ /= w;  *p++ /= w; *p++ /= w;
  355.     *p = 1.0;
  356.   }
  357. }
  358. #endif
  359.  
  360. /*
  361. ** Return the vector V=P1-P0.
  362. */
  363. extern VectorType4D Pdiff(P1,P0)
  364. PointType4D P1,P0;
  365. {
  366.   VectorType4D V;
  367.  
  368.   NormalizePoint(&P0);
  369.   NormalizePoint(&P1);
  370.  
  371.   V.V[0] = P1.P[0] - P0.P[0];
  372.   V.V[1] = P1.P[1] - P0.P[1];
  373.   V.V[2] = P1.P[2] - P0.P[2];
  374.   V.V[3] = 0.0;
  375.  
  376.   return V;
  377. }
  378.  
  379. /*
  380. ** Return the vector s*V.
  381. */
  382. VectorType4D Vscale(s,V)
  383. double s;
  384. VectorType4D V;
  385. {
  386.   V.V[X] *= s;  V.V[Y] *= s;  V.V[Z] *= s; V.V[W] = 0.0;
  387.  
  388.   return V;
  389. }
  390.  
  391. /*
  392. ** Return the vector V1+V1.
  393. */
  394. VectorType4D Vadd(V1, V2)
  395. VectorType4D V1, V2;
  396. {
  397.   VectorType4D V;
  398.  
  399.   V.V[X] = V1.V[X] + V2.V[X];
  400.   V.V[Y] = V1.V[Y] + V2.V[Y];
  401.   V.V[Z] = V1.V[Z] + V2.V[Z];
  402.   V.V[W] = 0.0;
  403.  
  404.   return V;
  405. }
  406.  
  407. /*
  408. ** Return the point P+V.
  409. */
  410. PointType4D PVadd(P,V)
  411. PointType4D P;
  412. VectorType4D V;
  413. {
  414.   PointType4D p;
  415.  
  416.   NormalizePoint(&P);
  417.   p.P[X] = P.P[X] + V.V[X];
  418.   p.P[Y] = P.P[Y] + V.V[Y];
  419.   p.P[Z] = P.P[Z] + V.V[Z];
  420.   p.P[W] = 1.0;
  421.  
  422.   return p;
  423. }
  424.  
  425. /* Return the negative of vector V */
  426. VectorType4D Vneg(V)
  427.   VectorType4D V;
  428. {
  429.     V.V[X] = -V.V[X];
  430.     V.V[Y] = -V.V[Y];
  431.     V.V[Z] = -V.V[Z];
  432.     /* since V is a vector, V.V[W] is always 0. */
  433.     return V;
  434. }
  435.  
  436. /* Returns the midpoint to P1 and P2 */
  437. PointType4D Midpoint(P1,P2)
  438.   PointType4D P1,P2;
  439. {
  440.       PointType4D M;
  441.  
  442.     NormalizePoint(&P1);  
  443.     NormalizePoint(&P2);
  444.     M.P[X] = (P1.P[X] + P2.P[X]) / 2.;
  445.     M.P[Y] = (P1.P[Y] + P2.P[Y]) / 2.;
  446.     M.P[Z] = (P1.P[Z] + P2.P[Z]) / 2.;
  447.     M.P[W] = 1.;
  448.     return M;
  449. }
  450.  
  451.